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■ ABSTRACT 

O 
SO 
O 

■ We present a detailed computational study of the assembly of protostellar 

On' 

^ ' disks and massive stars in molecular clouds with supersonic turbulence. We 

, follow the evolution of large scale filamentary structures in a cluster-forming 

si 

C3 ■ clump down to protostellar length scales by means of very highly resolved, 3D 

adaptive mesh refined (AMR) simulations, and show how accretion disks and 
' massive stars form in such environments. We find that an initially elongated 

cloud core which has a slight spin from oblique shocks collapses first to a 
filament and later develops a turbulent disk close to the center of the filament. 
The continued large scale flow that shocks with the filament maintains the high 
density and pressure within it. Material within the cooling filament undergoes 
gravitational collapse and an outside-in assembly of a massive protostar. Our 
simulations show that very high mass accretion rates of up to 10~ 2 Mq yr _1 
and high, supersonic, infall velocities result from such filamentary accretion. 
Accretion at these rates is higher by an order of magnitude than those found in 
semi-analytic studies, and can quench the radiation field of a growing massive 
young star. Our simulations include a comprehensive set of the important 
chemical and radiative processes such as cooling by molecular line emission, 
gas-dust interaction, and radiative diffusion in the optical thick regime, as well 
as H2 formation and dissociation. Therefore, we are able to probe, for the first 
time, the relevant physical phenomena on all scales from those characterizing 
the clump down to protostellar core. 
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1 INTRODUCTION 



Arguably the most comprehensive theoretical picture that we have of star formation is 
that stars are a natural conseque nce of supersonic turbulence within self-gravitating, molec 
ular clouds (see recent re views 



Mac Low fc Klessen 



Ballesteros-Paredes et al 



200 



El megreen fc Scalo 



2004; 



2004). It has been appreciated for some time that turbulence could 



play an important role in controlling th e formation of cosmic structure in general and stars 
in particular (e.g.. lvon Weizsackerlll951[ ). Supersonic turbulence is observed in most if not all 



giant molecular clouds (GMCs) and is important because it rapidly sweeps up large volumes 
of gas and compresses it into systems of dense filaments. The density jump that is associated 
with isothermal shock waves of Mach number M. scales as ~ A4 2 , so that this process may 
ulti mately be respo nsible for forming the dense clumps and cores seen in molecular clouds 



le.g.. 



Padoan 



1995). An attractive aspect of this theory is that since turbulence is associated 



with a broad spectrum of velocities, a large range of of physical scales and masses should 
arise for both the filamentary structure, as well as the properties of the dense clumps and 
cores that form within them. 

Observational surveys confirm that filamentary substructure characterizes the internal 
organization o f molecular clouds on man y scales. It is clearly evident in the Orion A molecu- 



lar cloud (e.g. 



Johnstone Rally 



19_9_9J) where, in addition to showing the obvious integral- 



filament shaped structure, one also sees sma . 



millimeter studies of Orion B (Mitchell et al 



ler st ructures of 1.3 pc in scale. JCMT sub- 



20011 ) reveal a plethora of filamentary struc- 



tures and their em bedded cores. This pattern is also seen in the observational studies by 



Fiege et al. 



(2004), wherein a few bright cores are seen to be embedded in a larger, isolated 
filamentary structure. Similar results are also seen in more embedded regions such as the 
Lupus 3 cloud where strong links between filaments and emerging star clusters are observed 
flTeixeira et aJbfloj) 

Computational studies of turbulence (e.g. 



Por ter et al 



1994) have made important 



progress over the last decade both in their dynamic range and in the list of physical processes 
that have been added to them. Simulations show that the shocks that are created in super- 
sonic turbulence produce a rich filamentary substructure. A broad spectrum of smaller, con- 
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densed regions is also produced (e.g., 



Klessen fc Burkert 



2000; 



Padoan fc Nordlund 



2002). 



Detailed comparisons with the observations show that there is a good correspondence be- 
tween numerically determined core mass spectra and the su bmillimeter surveys of star form- 
ing cores in molecular clouds (e.g. 



Tillev fc Pudritz 



20041 ) (henceforth TP04). One of the 



most important consequences is that the core mass spectrum (CMS) is found to be similar 
to that of the initial mass function (IMF) that characterizes the stellar mass spectrum . 
Cores in these numerical experiments are also born with a broad distribution of spins, which 
are a consequence of the oblique shocks that create them (e.g., TP04). The general velocity 
field created in simulations of compressible gas has also be en shown to haye a v ery good 



correspondence with that observed molecular cloud spectra (Fal garone et al 



1994- Finally, 



structure formation in supersonic flows occurs rapidly - typically in a sound crossing time 
or or two. This can in principal account for the rapid time-scales that characterize the for- 
mation of star clusters such as the Orion Nebula Cl uster (ONC) - who se members have a 
mean age < lMyr with a spread of less than 2Myr (iHillenbrandl 11997) . Thus, the univer- 
sal properties of turbulence potentially provides a comprehensive explanation for the broad 
range in stellar masses and spins that characterize star formation wherever we have been 
able to measure it. 

The turbulent fragmentation theory of star formation makes a much broader range of 
testable predictions than the older theory that featured the gravitational collapse of nearly 
spherical, isolated cores, often characterized as singular isothermal spheres. This latter pic- 
ture did not provide a clear physical explanation for the origin of the IMF. Moreover, the 
model predicted that accretion rates that are related to the cube of the (isothermal) sound 
speed c, i.e. M a ~ c 3 /G. While adequate to explain the formation of low mass stars in atypi- 
cal clouds such as Taurus, this accretion rate is orders of magnitude too small to account for 
the formation of h igh mass stars in ob served regions of clustered star formation such as the 



Orion Nebula (see 



Benther et al 



2006, for recent review). It is known that accretion rates of 



order 10 3 yr 1 can overcome the effects of radiation pressure and result in massive star 



formation (Wolfire fc Cassinelli 



1987). It has been suggested that the cores in which massive 



stars form can achieve these high rates essentially due to 



is found in such regions (e.g. McLaughlin fc Pudritz 



;he high turbu 



mi 



McKee fc Tan 



ent pressure that 



2003). 



While high pressure must play an important role in the formation of massive stars, these 
latter theoretical studies did not address how it is that low mass stars would form in the 
high pressure region that is producing the massive stars - both populations are co-spatial. 
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Moreover, it has been shown in simulations that the turbulence envisaged to drive a high 
mass accretio n rate onto a single massive star will instead fragment the gas into smaller pieces 



quite rapidly (jDobbs et al 



2005). Turbulence creates structure, and it is this structure, not 



pressure, that ultimately determines mass flow. 

A more fundamental shortcoming of the isolated core picture is that it fails to explain 
how stars can span more than 3 decades in mass. While the process of ambipolar diffusion - 
often invok ed in older mode ls as a regulator of gravitational collapse of isolated magnetized 



cores (e.g., 



Shu et al. 



19871 ) - certainly does mediate the eventual collapse of sufficiently 



massive cores, it does not provide a first-principals explanation for why such a wide range of 
stellar masses should exist in the first place. Similarly, Jeans mass arguments to "predict" 
stellar masses are not really very predictive when applied to media that are dominated by 
large density contrasts. As low and high mass stars form in clusters and are essentially 
co-spatial they are likely to acquire their masses by similar mechanisms. In both cases, we 
expect that filamentary accretion will play a role in the accumulation of the cores out of 
which they form. 

In this paper, we explore the multi-scale nature of the formation of filamentary large scale 
structure and how this controls the formation of protostellar disks and stars in turbulent 
cluster-forming clumps (these have typical sizes of < lpc). We accomplish this by using 



an Adaptive Mesh Refinement (AMR) code known as FLASH ((Frvxell et al 



2000). Our 3D 



hydrodynamical AMR simulations include many physical processes such as cooling due to 
dust, molecules, molecular hydrogen dissociation, as well as radiative diffusion out of opti- 
cally thick regions. We use as initial conditions numerical data from our recent simulations 
of cluster formation in molecular clumps (TP04) that show that a complete mass spectrum 
of gravitationally bound cores can be formed by supersonic turbulence within molecular 
clouds that closely resemble the core mass function of clouds such as Orion. 

We follow the filament that forms the first massive star in our simulation, and find 
that it has high density and dynamic pressure that is maintained by continued inflow and 
shocking into the filaments. The collapse of this dense material along the filament and into 
the disk and star is the key to understanding this problem. Lower mass stars in this picture 
also form in filamented structure, but which are less compressed and on smaller scales. Our 
simulations trace how the most massive star forming region is assembled and how collapse 
and the formation of protostellar disks can occur, by resolving the local Jeans length down 
to scales approaching that of the protostar itself. 
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We find that filaments play a dominant role in controlling the physics, accretion rate, 
and angular momentum of th e much smaller-scale accretion disk that forms within such 



Balsara et al 



200 lh . Large scale filamentary flows sustain 



collapsing structures (see also 
accretion rates that are orders of magnitude greater than both (i) the naive scalings that are 
derived from the virial theorem applied to uniform, 3D media, or even (ii) the collapse of 
isolated Bonner- Ebert spheres. The large scale structuring of molecular clouds into filaments 
therefore has profound effects on the rate of formation of disks and stars. 

The formation and subsequent fragmentation of large scale filaments can also lead to the 
formation of (ma ssive) binaries or multiple star systems which was shown in a ser ies of SPH 
simulations (e.g. 



Turner et al. 



1995; 



Whitworth et al. 



1995; 



Bhattal et al. 



19981) . In those 



studies, fragmentation is mainly the result of the complete formation of molecular hydrogen 
(which reduces the temperature in the cloud cores at high densities) during the gravitational 
collapse of the large scale filament. 

This paper proceeds by first describing our initial conditions and numerical procedures 
($2J. We then go on to systematically study the evolution of large scale filaments (fjHJ) that 
feed the formation of a massive star, the formation and evolution of an accretion disk (jSj), 
and finally the accretion of the central star (jjJJ) within the disk. We conclude our results in 
the discussion of (^BJ and give a compilation of our cooling approach in the appendix 1X1 



2 NUMERICAL METHOD AND INITIAL CONDITIONS 



Many studies of star formation have examined the self-similar collapse of highly idealized 
structures such as singular isothermal spheres. An initial state that is much closer to ob- 
servational reality is the 3D a Bonn er-Ebert sphere - who se c ollapse has received far les s 
theoretical attention until recently. In lBaneriee et al.l (|200J) and lBaneriee fc Pudritz (2006), 
we adapted the FLASH code for star formation applications and studied the three dimen 



siona l collapse of ro tating, magnetized and non-magnetized, Bonnor-Ebert-Spheres (jEbert 



12£5j 



Bonnor 



1956). We will compare our simulations with the results of the collapse of 



isolate B-E models later in the paper. As noted above however, isolated objects like this are 
not typical structures of observed molecular clouds. 

This paper attacks the broader problem of the physical connection between large scale 
structure formation (in this instance, on the scale of the cluster forming clump) and star 
formation. We use as initial conditions the final state of simulations of turbulent core for- 
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Figure 1. Initial setup as "seen" by the FLASH code. The low density (density shown in gray scale) areas are less refined due 
to the applied Jeans refinement criterion. The resolution is indicated by the square blocks where each 3D block has 8 3 grid 
points. 



mation presented by TP04. The goal of this latter study was to study the origin of the 
core mass function in cluster forming clumps, and to identify the physical conditions that 
best matched the observations. As documented below, we found that there was a particular 
set of conditions that produced realistic properties of molecular cloud cores and that also 
abounded in filamentary-like structures. 

T he TP04 simulations we re hydrodynamic and were performed with the ZEUS 3D 



code ( Stone Norman 



on a static grid and had to be stopped at the point when 



the density of any of formed cloud cores reached the critical density of the Truelove crite- 



rion (Truelove et al 
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19971 ). Using the FLASH code (IFrvxell et alihoooh wh ich is based on 



the PARAMESH adaptive mesh refinement (AMR) technique (jOlson et al 



1999f) allows us 



to follow the collapse of the TP04 cores over many orders of magnitude in density and length 
scale without under-resolving the physical relevant scales. 

In practice we took the raw simulation data (velocity and density) from one of the later 
output files in TP04, which are available in HDF format, and initialized the FLASH variables 
with these ZEUS-data. During the initialization we use a Jeans refinement criterion where 
we resolve the local Jeans length, 

„2 \ 1/2 



A. 



7T C 



Gp 



(1) 



by at least 8 grid points (c is the isothermal sound speed, p is the matter density, and 
G is Newton's constant). This density dependent initialization results in an initially non- 
homogeneous grid structure compared to the homogeneous ZEUS 3D grid. By de-resolving 
the low density regions in the simulations box we can save a substantial amount of compu- 
tational time without under-resolving the physical scales which are governed by the Jeans 
length during the collapsing stage. In Fig. ^ we show a 2D slice of the initial grid structure 
as "seen" by the FLASH code. 

For the subsequent simulation we increased the resolution of the Jeans length to at 
least 12 grid points. Note, that the Truelove criterion recommends that the Jeans length is 
resolved by only 4 grid points to avoid numerical fragmentation and ensure convergence. 



2.1 Initial physical parameters 

For the initial setup, we take the late stage numerical data of run B5 from the TP04 
simulations. This run has the following initial values for the side length of the simulation 
box, L = 9.8 x 10 17 cm = 0.32 pc, total mass, M tot = 105.1 Mq, isothermal sound speed, 
c = 0.408 km sec -1 , and a rms velocity, t> rms = 5 c with a Kolmogorov-type spectrum, respec- 
tively. The simulation run B5 resulted in a mass spectrum of bound cores that has a Salpeter 
slope at high masses (in agreement with the observations). At the very highest masses, one 
finds a few objects, the most massive of which was the first to collapse. This is the object 
whose collapse we continued to follow with the FLASH code 1 . At the time we re-started 



1 For better visualization we shifted the all simulation variables so that the peak density appears at the simulation center 
before we re-run the FLASH simulation. This could be done without changing any physical properties as both simulations use 
periodic boundary conditions. 
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the collapse simulation with the FLASH code the peak mass density is 1.35 x 10 -15 gcm -3 
which corresponds to an initial free fall time, t s - 



1 3tt I '32 Gp ~ 1800 yrs. The mass en- 
closed within a radius of 2 x 10 17 cm around the peak density, which is approximately the 
size of the cloud core, is ~ 23 Mq. This is also the Jeans mass of the initial TP04 simula- 
tion and represents a fairly massive cloud core mass and i ts ev olution can be comp a red to 
the collapse simulation presented in IBaneriee et al.l (J200J) and IBaneriee fc Pudrit z (2006) 
which are based on the phys ical properties of a 170 Mq cloud core and the 2.1 Mq cloud 
core Barnard 68 observed bv lAlves et al.l ((20011 1 . 



2.2 Cooling 



We extended the physical and chemical parameter space used in our earlier work to in- 
corporate the relevant cooling processes important for star formation in molecular clouds. 
These processes are summarized in complete detail in the appendix [X] to which we refer 
the reader for further information. We included cooling by molecular line emission, cooling 
from gas-dust interactions, radiative diffusion in the optical thick regime, and cooling by H2 
dissociation. 

For the remainder of this subsection, we give the more general reader an idea of how we 
proceeded (one may wish to go directly to the results in section^)). We first draw a distinction 
between the optically thin and optically thick regim es. For molecu l ar line emission in the 



optically thin 



Baneriee et al 



i mit, w e us ed the cooling rates d ue to 



optically thin regime was treated as in 



Goldsmith 



Neufeld et al 



(1995), as employed in 



(200 J) and iBanerjee fc Pudritz! (1200611 . Cooling by gas-dust transfer in the 



(|200l[ l. We show in the appendix that the 
dust temperature approaches the gas temperature at densities exceeding n > 10 7-9 under 
typical conditions. Even under these conditions however, we find that the dust is still a very 
efficient coolant. 

In the optically thick regime, the radiation from the dust surfaces (Eq. (|A4|1 ) cannot 
escape instantaneously and, by definition, the mean free path of the radiation field (~ l//t) 
becomes smaller than the typical core size of the dense core. For most collapse situations the 
core size is given by the core's Jeans length, Aj = yV c 2 s /Gn p coia - Therefore, we estimate 
the optical depth as 

t^kXj, (2) 
where k is the dust opacity. This approximation is in accordance with the results from 
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Silk 



opacit y limited fragmentation and star formation by lLow fc Lvnden-Belll (|1976f ) and 
( 19771 ) who use the expression above in their studies. 

Using the fact that the dust temperature gets close to the gas temperature (see section 
above) we can use Eqs. (jA5| ) and (0) to estimate the critical density at which the gas is not 
cooled efficiently, i.e. when r = 1 

rp x —5 



n r 



2.1 x 10 



10 



20 K 



cm 



(3) 



In 



Baneriee et al.l f)2004l ) we showed that this critical density can be used to estimate the 
distance from the center at which the first shock fronts build up. By updating our previous 
results with the above critical density we find 

r crit ~ 17 AU ( n ^ Q -Y 1/2 . (4) 

V2.1 x 10 lu cm~-V 

Test simulations of spherical collapse and the simulation reported here confirm the appear- 
ance of first shocks at distances close to the above value. 

We treated the propagation of radiation in the optically thick regime in the radiative 
diffusion limit (see the appendix). We also incorporated dust properties such as the frequency 
independence of cooling at temperatures above 100 K, as well as the melting of dust grains 
at temperatures around 1500 K. 

Finally, we included both the formation and dissociation of molecular hyd rogen. The for- 



matio n rate of H 2 molecules on dust surfaces was computed using the results of 



Hollenbach & 



McKef 



1979). We used the dissociation rates of molecular hydrogen as computed in lShapiro fc Kane 



1987). 



In order to understand the effects of these different processes in a dynamic calculation, we 
first ran a series of test simulations of the collapse of isolated, non-rotating, non-magnetized 
B-E spheres. Fig. |21 shows the result of our BE-collapse which include the process of H 2 dis- 
sociation (this simulation started with the assumption that the gas is completely molecular 
to begin with) . As shown in the appendix, the temperature in the radiation diffusion regime 
rises as T oc n 1 ^ 3 . At the density n cr ~ 10 16 cm~ 3 , the temperature reached ~ 1200 K, which 
is hot enough to effectively dissociate hydrogen molecules. However, the subsequent reduc- 
tion of the gas temperature reduces the dissociation efficiency leading to a self-regulated, 
slow, disintegration of hydrogen molecules over a wide range of densities. 

We now turn to an analysis of the results of our turbulence simulations and the formation 
of massive stars. 
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Figure 2. Shows the evolution of the hydrogen densities and gas temperature of our non-rotating, non-magnetized collapse 
simulation. As predicted from our theoretical considerations the dissociation of H2 molecules is a self-regulating process over 
a wide range of densities. The efficient cooling from the dissociation process (see Eq, IA28l ') reduces the gas temperature and 
subsequently slowing down the dissociation of H2 molecules. 

3 EVOLUTION OF THE LARGE SCALE FILAMENTS 

The initial state for our simulation, shown in Figure ^ (a 3D version of this is rendered in 
Figure 7 of TP04), already shows that large scale structure has formed in our simulation 
volume. Hundreds of lower mass cores are present in this volume, and as the virial analysis 
in that paper showed, many of the cores have started to collapse. The first core to violate 
the Truelove criterion in TP04 is the most massive core. As already noted, the oblique 
shocks which are responsible for the core formation are also responsible for the non-vanishing 
distribution of angular momenta in the cores. As we will see later, this initial spin results in 
a disk during the contraction phase. The AMR simulation zeros in on the evolution of the 
filamentary structure associated as it collapses to form a disk and protostar. The time-scale 
for the development of this nonlinear, gravitationally dominated, structure is very fast - 
amounting to several thousand years. In this time frame the larger filament scale remains 
fairly unaltered as material starts to drain into the assembling protostellar disk. 

The forming accretion disk is deeply embedded in the much larger scale filament as is 
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Figure 3. Density isosurfaces at different scales: 2400 AU, 150 AU, and 37 AU. The panels show the large scale filament, the 
protostellar disk within the filament, and the inner region of the disk. The density isosurface shown in opaque green in the 
left panel corresponds to the density shown in transparent green in the middle panel. Correspondingly, the opaque blue in the 
middle panel and the transparent blue in the right panel show the same density isosurface. 

nicely shown in Figure H3 The large scale filament is about 2600 AU in length and one sees 
the converging accretion flow on the 150 AU scale view presented in the middle panel. On 
this intermediate, 150 AU scale, the disk grows within a filament that is developing out of a 
sheet-like structure. As we will see, the transfer of material from the sheet into the filament is 
associated with net angular momentum - and tends to resemble the growth of a large-scale, 
spinning vortex. The disk, forming on yet smaller scales, acquires this angular momentum 
from this larger scale process. We see the actual disk at this time in the right panel, on a 
scale of only 37 AU in diameter. 

More details about the multi-scale structure of this filament are shown in the three 
Figures |U El These represent 3 different 2D cuts through the simulation data, at the end 
of our simulations. The structures are shown at increasingly resolved (finer) physical scales. 
Figure 0] and 03 feature cuts down the length of the filament which extend roughly along 
the vertical (z) axis ( xz and yz planes respectively). We clearly see a resolved filamentary 
structure on the largest scales (10 17 — 10 16 cm) scale in the upper left frame of the picture. 
On the smallest scales (highest resolution), shown in the bottom right panel, we see the 
formation of a disk. 

Two kinds of large scale velocity fields are shown in these figures: 

(i) flow of material into the filament from larger scales. This is a consequence of the large 
scale velocity field associated with the original supersonic turbulence field (most of the power 
is on the largest scales - as in any general turbulent flow); 

(ii) large scale flow of material along the filament that converges to form a disk near to its 
mass center. 
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Figure 4. Shows 2D slices of the matter density (logarithmic in grams per cubic centimeter) in the xz plane at different scales 
from the latest stage of our simulation. The panels a to / show the density structure on continuously smaller scales: 0.3 pc (a), 
8.2 X 10 3 AU (b), 1 X 10 3 AU (c), 128 AU (d), 16 AU (e), and 2 AU (f). These slices cut parallel through the filament and show 
the protostellar disk edge on. The arrows indicate the velocity field. 
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Figure 5. Shows 2D slices of the matter density (logarithmic in grams per cubic centimeter) in the yz plane at different scales 
from the latest stage of our simulation. The panels a to / show the density structure on continuously smaller scales: 0.3 pc (a), 
8.2 X 10 3 AU (b), 1 X 10 3 AU (c), 128 AU (d), 16 AU (e), and 2 AU (f). These slices cut parallel through the filament and show 
the protostellar disk edge on. The arrows indicate the velocity field. 
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Figure 6. Shows 2D slices of the matter density (logarithmic in grams per cubic centimeter) in the xy plane at different scales 
from the latest stage of our simulation. The panels a to / show the density structure on continuously smaller scales: 0.3 pc (a), 
8.2 X 10 3 AU (b), 1 X 10 3 AU (c), 128 AU (d), 16 AU (c), and 2AU (f). These slices cut perpendicular through the filament and 
parallel through the disk plane. The thin sheet attached to the filament is also clearly seen in the panels c and d. The arrows 
indicate the velocity field. 
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In the first case, the large scale flow continues to shock the filament and maintains the 
high pressure within it preventing its dispersal. We have measured the pressure in the initial 
state of this filament to be of the order of P/ks — 10 9 K cm -3 , which is high but expected in 
a converging flow at Mach numbers of order 5. In the second case, the material in the filament 
undergoes gravitational collapse along the filament and towards the point of concentration 
- the disk. 

We show the disk (edge-on) close-up in the bottom panels of Figs. 0] and |5J The filamen- 
tary collapse flow undergoes a shock with the dense and slowly cooling gas, and material 
joins the forming disk. Given that this is a highly asymmetric, 3D simulation, it is interest- 
ing that this small scale disk is not badly warped, although it cannot be described as being 
perfectly axisymmetric. 

A top-down view of the forming disk, and another aspect of the accretion flow that is 
giving its life, is seen in Figure El This cross-cut (xy plane) through the mid-plane of the 
disk shows that material from the large scale filament and the attached sheet is accreting 
onto the disk. Panel d of this figure especially shows an accreting structure that looks like 
a picture of accretion flow in a Roche-overflow in a binary system. The disk at the highest 
resolution (panel f) shows the presence of a bar, which is instrumental in driving a rather 
high accretion rate through the disk, as we shall see. 

While we do not follow the collapse of smaller regions in Figure Q in this paper, we know 
that filamentary structure is associated with low mass cores as well (TP04). These varied 
filaments are the result of shocks on different physical scales and different intensities. This 
implies that the mass flow rates through them and onto their embedded cores would be 
very different. Thus, low mass stars would also be expected to form through the filamentary 
accretion we study here. 

4 DISK FORMATION 

We present in Figures [7| and |S] a close-up snapshot of the temperature structure of the disk, 
on ~ 16 AU scales. Figure [7| shows the accretion shock very clearly. Material, raining into 
the disk from large scale flows along the filament, shocks with the disk because the gas 
cannot cool quickly enough at these densities. One notes that rather high temperatures of 
the order of 1000 K are associated with this accretion shock which envelops the entire disk. 
The post-shock gas cools rather quickly, and is found at much more modest temperatures 
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Figure 7. Temperature contours (in Kelvin) in the xz and yz plane which cut perpendicular to the disk plane (scale: ~ 16 AU). 
The first shocks appear in the region of the optical thin - optical thick transition. The arrows indicate the velocity field. The 
corresponding density maps are shown in panel e of Fig. HI and Fig. |S] respectively. 



of 200 - 400 K. As demonstrated in 



1995 



Tomisaka 



2002: 



Baneriee et al 



Matsumoto fc Tomisaka 



2004) and by others (e.g. 



Yorke et al 



2004) this cooler post-shocked gas gets re- 



accelerated and shocks again further into the gravitational potential resulting in a double 
(or multiple) shock structure of the accreting gas. 

A top-down close-up of the disk's temperature distribution is shown in Figure |H1 We 
clearly see a well resolved large scale spiral wave pattern with two spiral arms (on 16 AU 
scales), that appears to be attached to a bar at the center of the disk. This is the state 
of the disk at the earliest times, long before the Class phase has ended. The presence of 
this density wave is central to the rapid transport of disk angular momentum that we shall 
examine later. 

In order to understand what kind of emission one might expect from such a disk (in 
accordance with our radiation diffusion approximation), we present an optical depth map 
- which maps out values of the Rosseland mean opacity of the disk - is shown in the two 
panels of Figure OH We see the disk in edge-on as well as a top-down views. The post shock 
gas is very dense and the grains give it a high optical depth ranging from r ~ 10 3 — 10 6 in 
the dense regions on ~ 10 AU scales. One also sees that the spiral wave features seen in the 
temperature and density maps also show up in the optical depth maps of the disk's radial 
structure. 

The evolution of the radial structure of the disk is shown in Figures ITOl and HT1 2 . In 



2 These radial profiles show spherically averaged quantities with the center at the peak density 
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Figure 8. Temperature contours (in Kelvin) in the disk plane (xj/-plane) at ~ 16 AU scales. Clearly visible are the two arm 
spirals and the center bar in this temperature map. The corresponding density map is shown in panel e of Fig. [5] 



these very early and deeply embedded stages of disk and star formation, we find that the 
disk is actually more massive than the central protostellar object for a brief period. This 
is seen in the left panel of Figure wherein the mass inside of 1 AU at the latest time is 
roughly 10~ 2 Mq, whereas the mass enclosed within 10 AU is about 10 times larger; 10 _1 Mq. 
Another interesting feature of this mass distribution is that the mass within the disk and 
protostar is a small fraction of the total mass that is distributed on larger scales. At 10 17 
cm, this amounts to a total of 13 Mq while at a parsec scale, we have a total of 100 solar 
masses - reflecting the mass of the cluster forming clump. This distribution of mass is very 



Abel et al 



simila r in character to simulations of the first star formed in the universe, by 
(2002). Although the cooling function of gas forming the first star is - of course - vastly 
different than later epochs that we are simulating here, nevertheless the two results both 
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Figure 9. Optical depth (r logarithmic see Eq. J2J) perpendicular to the disk plane (left panel) and in the disk plane (right 
panel) at ~ 32 AU scales. The optical depth varies over six orders of magnitude within the protostellar disk and resembles its 
spiral structure. 

show that massive stars are like the tips of icebergs in much more massive and extensive 
mass distributions. 

The right hand panel of Figure El shows that the central density of the disk grows very 
rapidly with time. At each time, the run of density has a central, flat core-radius (which size 
is of the order of the core's Jeans length) with a steeply falling density dependence beyond. 
Variations in the power-law slope arise because of the onset of different cooling rates of 
material at different densities. 

The dominance of the disk over the protostellar mass at this early time does not violate 
any observations of disk-star systems since the latter are generally available typically for 
Class I TTSs and later. Theoretical calculations of disk formation often employ the assump- 
tion of self-similar collapse, in which the central object dominates the mass. We do not see 
this in our simulations. We also note that a dominant central mass does no t appear in our 



earlie r simulations of the hydrodynamic collapse of a Bonner-Ebert sphere (jBaneriee et al. 

2nd). 



Theoretical discussions of disk formation also often noted that disks needed to be lower 
in mass than their central stars, otherwise they would be out of equilibrium. This presents 
no problem to these early systems. In fact, rapid stellar formation would be highly favored 
if disks were more massive in their earlier stages. These highly unstable disks favor the 
appearance of high formation rates - which is what we will show later in this paper. 

We also do not see further fragmentation of the protostellar disk at this stage. The 
increasing pressure of this optical thick disk stabilizes the disk against fragmentation. This 
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Figure 10. Radial profiles of the enclosed mass and mass density at different times. The t = Oyrs graph shows the profile at 
the beginning of our simulation run and the subsequent labels indicate the time past since then. 
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Figure 11. Radial profiles of the surface density and temperature at different times. The t = Oyrs graph shows the profile at 
the beginning of our simulation run and the subsequent labels indicate the time past since then. 



shows again that the thermal evolution of the collapsing cloud core plays a crucial role in 
determining the fate of the protostellar system. For instance, molecular clouds will fragment 
easily in the regime wh ere abundant molecular hydrogen is formed and subsequently the 
pressure decreases (e.g., 



Turner et al. 



1995 



Whitworth et al. 



1995 



Bhattal et al. 



1998). 



Figure ITTI summarizes the evolution of the disk's radial column density and temperature 
distributions. The outer regions of the disk beyond 10 AU or so show a rather shallow column 
density profile, S oc r~ L25 compared to the minimum mass solar nebula. We also see that 
a high column density of ~ 10 3 g cm -3 is achieved at about an AU scale, which is then 
surpassed. We are charting the non-equilibrium early phase of disk formation and evolution 
just before the major delivery of mass into the central protostar has occurred. It is therefore 
not a surprise that transitory column densities in excess of the solar nebula model will arise. 
The radial temperature distribution in the disk scales as T oc r _1 at larger radii, with a flat 
distribution in the disk center. This steep temperature profile reflects again the early stage 
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during the assembly of the protostar where shocks from the infalling gas onto the adiabatic 
core are still present. It is interesting that while the envelop of this collapsing region features 
relatively smooth variations of the column and density profiles, the temperature profile falls 
steeply and then abruptly changes to the nearly isothermal background. This transition in 
disk temperature behavior may be the best way of detecting the outer accretion shock and 
hence the outer disk edge. 

The question of how angular momentum is distributed and transported in these early 
phases is arguably one of the most interesting and important issues for disk evolution, 
and massive star formation. We plot, in the left panel of Figure EE the specific angular 
momentum j(r) (angular momentum per unit mass) as a function of disk radius for different 
times. We see that for all times after the initial state, that j(r) is an increasing function 
of disk radius. Rayleigh's theorem informs us that this is in fact a stable distribution of 
angular momentum (rotating bodies with constant or increasing specific angular momenta 
are stable). Only one graph is different; this corresponds to the initial state where we see 
a leveling of j(r) at radii 10 15 — 10 16 cm. Since the MRI instability is not present here, we 
must conclude that the bulk of the disk must be relatively free from turbulence in this early 
stage. 

In the right hand panel of Figure ^1 we replot this data in another interesting way - now 
showing the specific angular momentum as a function of the mass enclosed. This is a good 
way of measuring h ow much angular momentum is extracted and transported away from the 



enclosed mass (e.g. lAbel et al.ll2002n . The figure shows that the specific angular momentum 



of material within the 0.01 Mq mass shell has quickly lost nearly 3 orders of magnitude of 
specific angular momentum, having fallen from an initial value of 10 20 cm 2 sec" 1 to about 
10 18 cm 2 sec -1 in only 2400 yr. This angular momentum transport is achieved by the spiral 
waves that we have seen in the earlier Figures (Figs. El El and El). As is well known, spiral 
wave torques are very efficient in transporting disk angular momentum because of their long 
lever arm. In the absence of either outflows or turbulence (through the MRI instability that 
would be active in the presence of a magnetic field), the plot shows that the non-equilibrium 
disk is very efficient in generating a spiral wave structure that efficiently transports angular 
momentum out of the central disk. 

We measure how far the disk is from a Keplerian distribution in Figure ED We normalize 
the measured rotation speed with the local Kepler velocity by using the mass in the interior 
radius. The left hand panel of this figure shows that this disk rotates more slowly by a factor 
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Figure 12. The evolution of the spherical averaged specific angular momentum as a function of radius (left panel) and enclosed 
mass (right panel). During the accretion phase angular momenta is exchanged between "mass shells". The t = Oyrs graph shows 
the graphs at the beginning of our simulation run and the subsequent labels indicate the time past since then. 
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Figure 13. The evolution of the rotational velocity compared to the Keplerian velocity as a function of radius (left panel) and 
enclosed mass (right panel). The t = Oyrs graph shows the profile at the beginning of our simulation run and the subsequent 
labels indicate the time past since then. 

of about 2 than Keplerian value within several 10s of AU. There are several peaks in the 
value of the rotation speed suggesting that shock waves within the disk are active and may 
be playing a role in helping to redistribute disk angular momentum. 

The right hand panel of Figure shows that the rotation speed of enclosed mass shells 
increases quickly over the simulated time period which for the 0.1 Mq mass shell doubles 
in spin from 1/3 to roughly 2/3 of its Kepler value. This reflects the growth in mass of the 
central protostellar core at the expense of disk material, over this time period. 



5 MASSIVE STAR FORMATION 

Although we stopped our simulation at a very early stage into the formation of the central 
star (pre-class object) we speculate that this object will grow quickly to form a massive 
star since it is embedded in a massive envelope and the mass accretion is sufficiently high to 
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Figure 14. The evolution of the radial infall velocity as a function of radius (left panel) and enclosed mass (right panel). The 
collapse proceeds from outside-in resembling a collapsing Bonner-Ebert sphere. The t = Oyrs graph shows the profile at the 
beginning of our simulation run and the subsequent labels indicate the time past since then. 



overcome the radiation pressure from the burning young massive star. Most of the material 
that is destined to form the massive star in this simulation is assembled by channelling the 
material flowing along the spinning filament, through the accretion shock and into the star 
by accreting through the disk. Of greatest importance is the question of just how large the 
disk accretion rate is in these early phases. 

We plot the spherically averaged, radial infall velocity as both a function of radius, as 
well as of enclosed mass, in Figure El (left and right panels respectively) . The left hand 
panel shows some of the characteristics of outside-in collapse that distinguish the collapse of 
Bonner-Ebert spheres. In the left hand panel of this Figure, we see that the core is already 
in a state of infall in the initial state of our simulation. The outer envelope has begun to 
collapse while most of the interior material is still nearly at rest. As the collapse develops, 
the core radius progressively shrinks as its density grows. At the same time, the maximum 
infall velocity continues to grow. The outside-in velocity structure is essentially the same as 
that seen in the collapse of Bonner-Ebert spheres, wherein the ever increasing density of the 
core is supplied by collapsing material in the envelop. 

The collapsing gas eventually reaches a density, on small enough scales, that cooling 
is too slow and a shock is formed as the outer material arrives at the position of the more 
slowly contracting denser material. This happens first at a radial scale of a few 10 14 cm. This 
is th e well known appearance of the first core, whose interior contracts almost adiabatically 



ve-g- 



Larson 



1969). The 2D spatial cuts shown in Figure U\ clearly show that this is indeed 



the scale of the first shock, where dust cooling in the optical thick regime has become too 
inefficient to prevent the appearance of a shock wave. This change in the cooling behavior 
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Figure 15. The evolution of the mass accretion through spherical shells around the core center. The magnitude of the mass 
accretion is shown in absolute values and scaled to the isothermal quantity c? so /G (left panel), and scaled to the local cf ocgl /G 
accretion measure (right panel). The t = Oyrs graph shows the profile at the beginning of our simulation run and the subsequent 
labels indicate the time past since then. 



of the gas also appears as the change in density behavior seen in the left panel of Figure fTTl 



This scale can be predicted theoretically from 



m previous papers (see e.g. 



Baneriee fc Pudritz 



re co oling function itself, as we have noted 



2006|) 



Another interesting aspect of the collapse is that the infall speed is much higher than 
the original isothermal sound speed that characterizes the initial state, and reaches values of 
Mach 6 on 10 15 cm scales. It is also higher than the asymptotic infall velocity of collapsing, 
isothermal Bonner-Ebert sphere where the Mach number reaches ~ 3. The reason for the 
exceedingly high infall rate in our case is the initial supersonic turbulence. The right hand 
panel of Figure El shows the infall speed as a function of the mass enclosed. Only a small 
fraction of the infalling material reaches such high Mach number, typically less than 10% at 
these early times. 

Perhaps the most important aspect of filamentary accretion is in driving high accretion 
rates of material onto the disk. We show the spherically averaged accretion rates as a func- 
tion of disk radius in the left panel of Figure El i n both absolute terms, and then scaled 
with respect to the initial isothermal sound speed, c iso , of the gas at the beginning of the 
simulation. As expected from the velocity behavior of the collapse shown in the previous 
Figure El the accretion rate peaks at some radius, and this peak moves inwards with time 
to smaller and smaller scales. The initial state is already collapsing, as already noted and 
expected from the TP04 simulation. The peak of the mass accretion rate is located at the 
same radius at which the infall speed achieves its highest value (compare left hand panels of 
Figs. 1141 andlT3|). Note that at 1 AU the infall rate is several lO _3 M0yr _1 which in the last 
frame is 10 times higher on sub AU scales. In terms of the initial isothermal sound speed 
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of the simulation, the infall accretion rate has reached more than 10 3 c 3 /G! This val ue is 



abou t 1000 times higher than predicted by the self-similar solution of singular collapse (Shu 



19771 ) and ab out 20 times hi gher than on e woul d expect for Bonner-Eb ert collapse, as first 



Larson 



9 ) 



(1969) see also 



Hunter 



1977J). 



calculated bv IPenstonl (J1969J) and 

The right hand panel of Figure Hoi once again plots the accretion rate as a function of disk 
radius, except normalized with respect to the local sound speed, c 3 cal /G. This latter quantity 
is certainly not isothermal, and is the sound speed in the local gas that is undergoing the 
many cooling processes that we discuss in the appendix. This local sound speed is of course 
greater than the initial isothermal sound speed that we start the simulation with because the 
gas heats substantially as a consequence of shocks. Thus, as the temperature rises towards 
the disk's interior the local sound speed increases as does the value of cf ocal /G. In the last 
time plotted, the ratio of the infall rate to that computed from this local sound speed is only 
a factor of 10 larger at scales around one AU. 



6 DISCUSSION AND CONCLUSIONS 

Our high resolution simulations of a cluster forming molecular clumps show that the forma- 
tion of large scale filaments plays a major role in the formation of accretion disks and massive 
stars. Our AMR technique tracked the collapse of the most massive structure which was the 
first to form in the cluster-forming clump which is our initial state. The rapid assembly of 
a massive star - at accretion rates that are an order of magnitude greater than analytic 
estimates - arises for 2 basic reasons: (i) supersonic turbulence rapidly sweeps up dense, 
large-scale structures by compressing a large volume of gas into sheets and filaments, and 
(ii) the high accretion rates in filaments that features the flow of some material at unusually 
high Mach numbers. 

6.1 Collapse in filaments 

The cooling function for gas in the initial state dictates that the dusty molecular gas remains 
essentially isothermal until densities of the order 10 10 — 10 11 cm -3 . The shock waves enor- 
mously compress the material, by factors of order M 2 where Ai is the Mach number. For 
our simulation, the initial rms Mach number (characterizing the entire turbulent spectrum) 
is 5, so that the forming molecular core in this filament has already had an enormous boost 
in its background density, and therefore background pressure. 
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The increase in gas density brought about by the shock then increases the infall rate in 
the filament by a large factor. Since the infall speed scales as v g oc p 1 ^ 2 , we anticipate that 
gravitational collapse in such a filament is boosted by a factor of v B oc M. over typical "virial 
theorem" estimates. 

We have also observed that the infall occurs in a fashion somewhat like the outside-in 
collapse of a Bonner- Ebert sphere. It is important to remember that a B-E sphere collapse has 
an infall rate that is much larger than given only by the sound speed; one finds M ~ / c 3 /G 
where / ~ 20 — 50 as was restated by Hunter (|l977 ). We argue that the factor should scale 
as / ~ A4 3 , i.e., that one replaces the sound speed in the virial formula with the free-fall 
speed. Thus, since the detailed numerics show that a maximum infall speed of M. ~ 6 is 
reached, then / ~ 216. This value gets us into the right regime but is still not the factor of 
10 3 enhancement that we see in the most extreme case in our. The difference could be due 
to a combination of factors including a geometric factor. 

The point is that the filamentary collapse delivers an accretion rate that is greatly in 
excess of what is needed to create a massive star in the conditions, on time scales of only 
a few thousand years from the point where these filamentary initial conditions have been 
created. 



6.2 Role of cooling 

Cooling plays an important role in these processes, because it ultimately controls how dense 
structures can become. As long as the gas is isothermal, it can become highly compressed as a 
consequence of shock waves. This early isothermal phase lasts until the gas starts to become 
optically thick to emission from dust, at densities exceeding 10 11 cm -3 . At this density, gas 
becomes optically thick to dust cooling and begins to contract quasi-adiabatically. This scale 
is characteristic of the first shock that appears at 10 — 20 AU (see Figure |7|) that defines the 
disk envelope. 

The second density of importance is at 10 16 cm -3 at which molecular hydrogen begins 
to dissociate. We find that an interesting, self-regulated state arises in which the gas is 
kept just cool enough to avoid both the total dissociation of molecular hydrogen, as well as 
the preservation of dust grains from evaporating. Molecular hydrogen cooling regulates the 
gas temperature. The dissociation of molecular hydrogen is i mportant bec ause it allows the 



collapse to resume, as has been noted by many authors (e.g. lLarsonlll969l ). 
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6.3 Conclusions 

Our basic conclusions are: 

• The assembly of protodisks and protostars is a rapid process in a supersonic, turbulent 
environment. 

• The most massive object forms first within the molecular cloud was the first to violate 
the Truelove condition in our simulation. This raises the question of how feedback from this 
(massive) star will affect the collapse and star formation in the many other smaller cores 
that formed in our simulation. 

• The collapse proceeds from outside-in where the accretion through the large scale struc- 
ture (a filament in our case) "feeds" the dense protostellar structure within it. 

• The initial spin of the filament from oblique shocks lead to disk formation because of 
efficient dust cooling until densities of ~ 10 11 cm" 3 . 

• Even in the absence of magnetic fields and/or outflows angular momentum is transfered 
very efficiently through spiral waves in the disk. 

• We observe very high accretion rates (M ~ 1CT 2 Mq yr _1 ) due to filamentary accretion 
of the supersonic gas onto the forming disks. These rates are 10 3 times larger than predicted 
by the collapse of singular isothermal spheres. We find that a reasonable scaling for filamen- 
tary accretion is M ~ fc 3 /G where the prefactor should scale as / ~ At 3 , i.e., that one 
replaces the sound speed in the virial formula with the free-fall speed. Even this does not 
quite account for the full effect, and there is an additional factor of order a few that likely 
is geometric in nature and accounts for filamentary, rather than spherical infall geometry. 
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APPENDIX A: APPENDIX: COOLING 

In this appendix we summarize the cooling processes which we included in our simulation. 



Al Molecular cooling 



As in 



Baneriee et al. 



(2004|) and 



line emissions calculated by 



Baneriee 



Neufeld et a 



fc Pudrita (|2006| ) we include cooling by molecular 



( 199 51) - These proces ses are most important in the 



Baneriee et al.l (J2004J) we showed that the collapse 



low density/low temperature regime. In 
will proceed isothermally until n ~ 10 7 5 cm -3 if only molecular cooling is considered. 



A2 Dust cooling 

Co oling by gas-dust energy transfer in the optical thin regime. Here, we follow the treatment 



by 



A 



Goldsmith! 1)20011 ) with the gas-dust energy exchange rate 

2 , A r~r~i - . m .1/2 



2 x 10 



-33 



/ AT\ / T. 



ergs cm 3 sec 1 



(Al) 



V cm- 3 J V K J U0K, 
where AT = T gas — T dust . This cooling process is very efficient until the core becomes optically 
thick. 

We calculate the the dust equilibrium temperature at each time step by solving 



I^cr -)- A gd A. 







(A2) 



to get the self-consistent dust temperature. In the above equation r cr and A dust are the 
heating by cosmic rays and the energy loss by dust (black-body) radiation, respectively. As 
long as AT > the dust is als o heated by the gas-dust interactions. The heating rate by 
Goldsmith! ( 2001 ) 



cosmic rays is 



3.9 x 10 



-28 



n(H 2 



cm 



ergs cm 3 sec 1 



(A3) 



where we assume a shielding parameter of \ = 10 4 . 
A2.1 Optical thin regime 

The dust grains loose their thermal energy by black-body like radiation: 



A ri „ s t — KG T,, 



(A4) 
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where k is the opacity (in cm -1 ) and a — f a is, up to a facto r of order unity, the Stefan- 



Boltzmann constant a. Using the values from iGoldsmithl (2001) 



k = 3.3x10- (^4) fJ^Vcm- 1 (A5) 



v cm- 3 J V 18.24 K, 

a = 6.85 x lO^ergscnT 2 sec" 1 K" 4 (A6) 
= 1.209 a 



the dust energy loss is 



A dust = 6.8 x 10- 33 (^j ergs cm" 3 sec" 1 (A7) 

Solving the dust equilibrium equation (|A2|) we find that the dust temperature, T dust , 
approaches the gas temperature, T gas at densities n(H 2 ) ^ 10 7-8 cm -3 depending on the initial 
gas temperature. 

When the dust temperature is close to equal to the gas temperature the effective gas 
cooling becomes 

¥ A dust (T dust = T sas ) (A8) 

where the dust temperature is replaced by the gas temperature. In this regime on can 
estimate the effective equation of state (EoS) by equating the compressional heating with 
the above cooling 

= ^t (A9) 



J- gas 
t-fl 



'if 



where t B — J 3%/ 32 Gn p is the free fall time and T dust is replaced by T gas in n. From Eq. ()A9jl 
we see that n oc T 10 which results in an effective adiabatic index, i.e. 7 cft = 1 + dlnT/dlnn, 
of 

Teff = ^ WneI1 T dust ~ T gas . (A10) 

This result shows that the dust is still a very efficient coolant even if the dust temperature 
becomes almost equal to the gas temperature. We summarize the temperature trajectory 
and the evolution of the effective equation of state in the Figures ED and IA1I which show 
the result of one of our non-rotating, non-magnetized, spherical collapse simulations. The 
simulation results confirm the existence of the transition region from the isothermal collapse 
to 7 cff = 1.1 at densities n > 10 9 cm -3 . 
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Figure Al. Shows the effective EOS as a function of density from our non-rotating, non-magnetized spherical, collapse 
simulation in the case where dust cooling and H2 dissociation is included. Clearly, two distinct transitions can be identified: 
at a density of ~ 10 11 cm -3 the core becomes optically thick and contracts on quasi-static trajectory, and at ~ 10 1 
dissociation becomes effective whereby slightly cooling the core region. 



A2.2 Optical thick regime 

We treat the propagation of the radiation field in the optical thick regime, i.e. r > 1 (see 
Sec. 12.21 for our definition of the optical depth r), as a diffusion process using the following 
approximation. The energy loss carried by the radiation flux F is 



A, 



VF. 



(All) 



In the Eddington approximation the flux can be written as fe.g.. iTurner fc Stonel l2001) 



F = -—VE 

3 K 



(A12) 



were c is the speed of light and E = a rad T 



la 



dust V^rad 



7.565 x 10 15 ergs cm 3 K 4 ) is 



the radiation energy associated with the dust. 
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In accordance with the approximation of the optical depth, Eq. (J2J), we use the Jeans 
length, Aj, as the typical length scale over which the the radiation field changes significantly. 
This allows us to approximate V ~ f/Xj and write Eq. (|A11J) as 3 
4 <tT 4 

A rad « - f — (optical thick regime) . (A13) 
3 k A j 

We note that the local treatment of the radiation energy loss, Eq. (jA13|) . is a good approxi- 
mation in an ongoing collapse situation where the Jeans length is the dominant scale of the 
system. The approximation becomes less accurate in a system which is governed by different 
independent length scales, like an equilibrium disk with a scale height h. 

The expression in Eq. f|A13|) sets an upper limit of the radiation loss in the optical thick 
regime and the dust cooling (Eq. ()A4jl ) is determined by 

A dust = A rad (A14) 

which leads to the replacement of the dust opacity in the optical thick regime 

k — > f I A j (optical thick regime) . (A15) 



where / = ^4cr/3a / « 1.05 /. 

In summary, we determine the dust temperature in the optical thin and optical thick 
regime by solving Eq. (jA2|) (using an iterative Newton method) with 

Adust = K cff & T dust (A16) 
with 

k cS = min (k, f Xj 1 ^ . (A17) 

The cooling in the optical thick regime is less efficient than in the optical thin regime. 

Therefore, the timescale which governs the evolution of the core becomes larger than the 

dynamical time scale of the core region and is given by the cooling time 

I n kT , k , 

- + - 2 (A18) 



aT*f/Xj 
n(H 

10 10 cm- 3 / VlOOK 



~ , / n(H 2 ) \ 1/2 / T \" 5 / 2 



3 The factor / could be a function of density and temperature which can be determined by comparing simulation results 
using the diffusion approach, i.e. Eqs. IA11I and IA12I . with results using Eq. JaTH . Here, we use a constant value which 
we determine by the comparison of Eq. IA18I with the results of our collapse simulation of a non-rotating, non-magnetized 
Bonnor-Ebert-Sphere. 
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In Fig. IA2I we show the density dependence of the core evolution time, t ovo i, sim = n/n 
from our simulation and the timescale from Eq. (jA18j) . Both timescales have the same 
dependence on density (oc n" 1 ^ 3 ) in the optical thick regime and the same amplitude if we 
choose / = 1.67. 

Equation (jA18|) does not constrain the density-temperature relation (i.e. the equation of 
state). But as the system collapses on a trajectory at which the thermal pressure is (almost) 
in balances with gravity this trajectory will be the path at which 

7e ff = Tent = 4/3 . (A19) 

Figs. IA1I and 121 show the evolution of the effective EoS (7 cff = dlnp/dlnp) and the 
temperature of the core region of our spherical collapse simulation. The collapse in the 
optical thick regime follows the path where 7 cff ~ 7 crit . 
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A3 Warm temperature regime and dust melting 

At temperatures of about 100 K the dust opacity becomes almost independent of the fre- 
quency and at temperatures around 1500 K the dust opacity dec reases sharply due to the 



Semenov et al 



2003). To 



melting of dust grains (for a recent compilation of opacities see 
incorporate these effects we fit the dust opacity by a piecewise power law (see also Eq. (|A5|) ) 

T 2 ; T < 200 K 

200K<T< 1500K (A20) 
T > 1500 K 



Jl-12 



A4 Ho formation and dissociation 



We use the iHollenbach fc McKeel (jl979f ) rates to account for the formation of H2. These 
authors compute the formation of hydrogen molecule on dust grain surfaces from collisions 
with atomic hydrogen. The formation coefficient, Rf, can be written as: 

T l/2 f 

-18 1 J a 



f 3X10 1 + 0.04 (T + T dust y/ 2 + 2 x 10~ 3 T + 8 x 10~ 8 T 2 
where T and T duat are given in K and the factor /„ is given by 

f 1 

J a 



(A21) 



(A22) 



1 + 10 4 exp (-600/T) ' 
To acco unt for the proces s of H? dissociation in the hot dens region we use the dissociation 



rates from 



Shapiro 



fc Kangi ()1987f ) which are recalculated functions based on the work by 



Lepp fc Shulll ()l983h . These dissociation rates are calculated from the collisional breakup of 

molecular hydrogen out of all 15 vibrational levels. 

The dissociation rate coefficient can be well fitted to the form 
log k H /k L 



log k D (n,T) = \ogk H 



(A23) 



1 + n/n ir ' 

where ko is given in units cm 3 sec -1 and ku and k^ are the coefficients for the low and high 
density regimes, respectively 

k L = 1.12 x 10~ 10 exp (-7.035 x 10 4 K/T) cm 3 sec" 1 ; H - H 2 (A24) 

= 1.18 x 10~ 10 exp (-6.950 x 10 4 K/T) cm 3 sec -1 ; H 2 - H 2 
k H = 1.20 x 10~ 9 exp (-5.24 x 10 4 K/T) cm 3 sec -1 ; H - H 2 (A25) 

= 1.30 x 10~ 9 exp (-5.33 x 10 4 K/T) cm 3 sec" 1 ; H 2 - H 2 . 

The transitional density n tI in units of cm -3 is given by 

logn tr = 4.00 - 0.416 x- 0.327 x 2 ; H - H 2 (A26) 
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= 4.845 - 1.3 x + 1.62 x 1 

with 

x = logT/10 4 K. 



Hr, - Ho 



(A27) 



The binding energy of H 2 molecules is 4.48 eV. For each collision with molecular hydrogen 
which leads to its dissociation this energy is lost from the (hydrogen) gas. Therefore, the 
effective cooling rate due to hydrogen dissociation is given by 



A diaa = 4.48eVnn(H 2 ) k D 



(A28) 



where n is either n(H 2 ) or n(H) depending on the collisional partner. 

The FLASH code provides data arrays which store the mass, X iy fraction of different 
species. We use this ability to keep track of the number densities of hydrogen in its different 
forms (i.e. molecular or atomic). The number density of the species % is given by 

= Xt , (A29) 

fii m H 

where /i, is the atomic weight of the species i, p is the total matter density, and m# is the 
hydrogen mass. 

We point out that we did not include effects from H 2 formation by three body processes 
which could still be important in the high density (> 10 16 cm -3 ) and warm temperature 
(< 1100 K) regime. We might therefore slightly overestimate the cooling ability from H 2 
dissociation. Three body for mation rate s are rather uncertain an d different compilations 



can be found, for instance, in 



Abel et al 



( 20021: IPalla et al.1 (jl983h. 
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